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Abstract. The method of moments in the context of Nonlinear Schrodinger Equations relies on defining 
a set of integral quantities, which characterize the solution of this partial differential equation and whose 
evolution can be obtained from a set of ordinary differential equations. In this paper we find all cases in which 
the method of moments leads to closed evolution equations, thus extending and unifying previous works in the 
field of applications. For some cases in which the method fails to provide rigorous information we also develop 
approximate methods based on it, which allow to get some approximate information on the dynamics of the 
solutions of the Nonlinear Schrodinger equation. 
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1. Introduction. Nonlinear Schrodinger (NLS) equations appear in a great array of con- 
texts PP as for example in semiconductor electronics |21 E], optics in nonlinear media 0], 
photonics [S], plasmas fundamentation of quantum mechanics [7], dynamics of accelerators 
mean- field theory of Bose-Einstein condensates |H] or in biomolecule dynamics jJOI- In 
some of these fields and many others, the NLS equation appears as an asymptotic limit for a 
slowly varying dispersive wave envelope propagating in a nonlinear medium 

The study of these equations has served as the catalyzer of the development of new ideas or 
even mathematical concepts such as solitons ^2 or singularities in EDPs I14| . 
One of the most general ways to express a nonlinear Schrodinger equation, is 

1 

i— = --/:^u + V{x,t)u + g{\u\^ ,t)u - i(T{t)u (1.1) 

where A — d"^ /dx\ + • • • + d"^ jdx\ and u is a complex function which describes some physical 
wave. We shall consider here the solution of on M" and therefore u : R" x [0, T] C, with 
initial values u{x,0) — uo{x) € X, X being an appropriate functional space ensuring finiteness 
of certain integral quantities to be defined later. 

The family of nonlinear Schrodinger equations contains many particular cases, depend- 
ing on the specific choices of the nonlinear terms g{\u\'^,t), the potentials V{x,t), the dissipa- 
tion (7{t) and the dimension of the space n. The best known cases are those of power-type 
gdup) = a\u\P or polinomial gdup) = ai|u|^^ -|-a2|u|^'^ nonlinearities. 

The potential term V{x, t) models the action of an external force acting upon the system and 
may have many forms. Finally, we include in Eq. Hl.l(l a simple loss term arising in different 
applications In many cases these losses are negligible, i. e. cr = 0. 
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The description of the dynamics of initial data ruled by this equation is of great interest for 
appHcations. Nevertheless, gathering information about the solutions of a partial differential 
equation that is nonlinear like constitutes a problem that is a priori nearly unapproach- 
able. For this reason, most studies about the dynamics of this type of equation are exclusively 
numerical. The rigorous studies carried out to date concentrate on (i) properties of stationary 
solutions 52], (ii) particular results on the existence of solutions and (iii) asymptotic 

properties P^ITl] . 

Only when n — l,g{\u\'^,t) — \u\'^,V{x,t) = 0, cr = it is possible to arrive at a solution of 
the initial value problem by using the inverse scattering transform method |12| . In this paper 
we develop the so-called method of moments, which tries to provide qualitative information 
about the behavior of the solutions of nonlinear Schrodinger equations. Instead of tackling the 
Cauchy problem Hl.l|) directly, the method studies the evolution of several integral quantities 
(the so-called moments) of the solution u{x,t). In some cases the method allows to reduce the 
problem to the study of systems of coupled ordinary nonlinear differential equations. In other 
cases the method provides a foundation for making approximations in a more systematic (and 
simpler) way than other procedures used in physics, such as those involving finite-dimensional 
reductions of the original problem: namely the averaged Lagrangian, collective coordinates 
or variational methods |2()l I21| . In any case the method of moments provides an information 
which is very useful for the applied scientist, who is usually interested in obtaining as much 
information as possible characterizing the dynamics of the solutions of the problem. 

It seems that the first application of the method of moments was performed by Talanov P2) 
in order to find a formal condition of sufficiency for the blowup of solutions of the nonlinear 
Schrodinger equation with ^dwp) = — [ttp and n — 2. Since then, the method has been applied 
to different particular cases (mainly solutions of radial symmetry in two spatial dimensions), 
especially in the context of optics where many equations of NLS type arise j23| . 

In previous researchs, the method of moments has been studied in a range of specific situ- 
ations but in all the cases the success of the method is unrelated to a more general study. In 
this paper we try to consider the method systematically and solve a number of open questions: 
(i) to find the most general type of nonlinear term and potentials in Eq. for which the 

method of moments allows to get conclusions and (ii) to develop approximate methods based 
on it for situations in which the moment equations do not allow to obtain exact results. 

2. Preliminary considerations. Let us define the functional space Q{H) as the space 
of functions for which the so called energy functional 

E{u)^{u,Hu)mu^^+ f G{\u\^t)dx (2.1) 

is finite, G being a function such that 9G(|up, t)/9|wp = g{\u\'^,t), {■,-)l^(R") denoting the 
usual scalar product in L^(R") and H — — -I- V{x,t). For the case V{x,t) = and g 
independent of time, several results on existence and unicity were given by Ginibre and Velo 

CHI 

As regards the case V ^ 0, which is the one we are mainly interested in our work, the best 
documented case in the literature is that of potentials with bounded in R" for every 

|a| > 2; that is, potentials with at most quadratic growth. Y. Oh proved the local existence 
of solutions in L^(R") and in Q{H) for nonlinearities of the type g{u) — —\u\p,0 < p < 4/n. 
However the procedure used allows to substitute this nonlinearity with other more general 
ones. It seems also posible to extend the results to the case in which the potential depends on 
t. 

Therefore, from now on we shall suppose that the nonlinear term satisfies the conditions set 
by Ginibre and Velo and that jD^yj is bounded in space for all \a\ > 2. Under these conditions 
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it is posible to guarantee at least local existence of solutions of Eq. Hl.l|l in appropriate 
functional spaces. 

2.1. Formal elimination of the loss term. In the first place we carry out the trans- 
formation |il5| 



u{x,t) = u{x,t)e'^o'''-^)d^, 



(2.2) 



which is well defined for any bounded function a{t) (this includes all known realistic cases 
arising in the applications). The equation satisfied by u{x,t) is obtained from the following 
direct calculation: 



at ot 



1 



Am + V{x, t)u + g{\u\^,t)u - ia{t)u 



= --Au + V{x,t)u + g{\u\'^,t)u 



where 



(2.3) 



From here on we shall consider, with no loss of generality, that a{t) = in Eq. assuming 
that this choice might add an extra time-dependence to the nonlinear term. 

3. The method of moments: GeneraUties. 

3.1. Definition of the moments. Let us define the following quantities: 

2 



h,jit) = / x''\u{x,t)\^dx 



k/2 
Xj U 



du{x,t) _du{x,t) 
dxj dxj 



dx. 



du{x, t) 



dxi 



dx 



Jit) = I Gi\uix,t)\'',t)dx, 



du 



dxi 



(3.1a) 
(3.1b) 
(3.1c) 
(3.1d) 



with j = l,...,7i and k ~ 0,1,2,..., which we will call moments of u{x,t) in analogy with 
the moments of a distribution. From now on and also in Eqs. (|3.1|) it is understood that all 
integrals and norms refer to the spatial variables x E M" unless otherwise stated. In Eqs. (|3.1|l 
we denote by u the complex conjugate of u. 

In some cases we will make specific reference to which solution u of Eq. (|1.1|) is used to 
calculate the moments by means of the notation: , etc. 

The moments are quantities that have to do with intuitive properties of the solution u(x, t). 
For example, the moment /o.o is the squared L^(R")-norm of the solution and therefore mea- 
sures the magnitude^ quantity or mass thereof. Depending on the particular context of appli- 
cation, this moment is denominated mass, charge, intensity, energy, number of particles, etc. 
The moments Iij{t) are the coordinates of the center of the distribution u, giving us an idea 
of the overall position thereof. The quantities l2.j are related with the width of the distribution 
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defined as = {J^„{xj — Ii^j)^\u\'^dx)^/^ ~ {h^j + I^jIqa — ^if j)^/'^ , which is also a quantity 
with an evident meaning. 

The evolution of the moments is determined by that of the function u{x,t). From now on 
we will assume that the initial datum uq (x) and the properties of the equation guarantee that 
the moments are well defined for all time. 

3.2. First conservation law. It is easy to prove formally that the moment /q.o is in- 
variant during the temporal evolution by just calculating 

^uAu - V{x, - dx 

-^uAu + V{x,t)\uf + g{\uf,t)\uf^ dx 

— (uAu — uAu) dx 
2 ^ ^ 

"^(/ iVMpdx- /" \\/u\^dx)=0, (3.2) 

where we have performed integration by parts and used that the function u and its derivatives 
vanish at infinity. 

Obviously the above demonstration is formal in the sense that a regularity which we do 
not know for certain has been used for u. Nevertheless, this type of proof can be formalized 
by making a convolution of the function u with a regularizing function. The details of these 
methodologies can be seen in ^S] or In this paper we will limit ourselves to formal 

calculations. 

4. General results for harmonic potentials. 

4.1. Introduction. From this point onward we will focus on the particular case of interest 
for this study when V{x,t) is a harmonic potential of the type V{x,t) — ^{x, A{t)x), where 
A is a real matrix of the form Aij{t) = X'f{t)Sij, with Ai > for i = 1, ...,n and Sij is the 
Kronecker delta. Bearing in mind the results of Section 12.11 the NLS equation under study is 
then 

" I {j^^ ^^'"^^^ j " ^^-^^ 

This equation appears in a wide variety of applications such as propagation of waves opti- 
cal transmission lines with online modulators [261 1271 128L EH] , propagation of light beams in 
nonlinear media with a gradient of the refraction index |3U[ I31| , or dynamics of Bose-Einstein 
condensates [Hj- Generically it can provide a model for studying some properties of the so- 
lutions of nonlinear Schrodinger equations localized near a minimum of a general potential 
Vix). 

4.2. First moment equations. If we differentiate the definitions of the moments Iij 
and Vo J , we obtain, after some calculations, the evolution equations 



dt 



= Vo,j, (4.2a) 



A?/i,„ (4.2b) 
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so that hjyj — 1, ...,n satisfy 







(4.3) 



with inital data /ij(0), /i.j(O) = Vb.j(O). These expressions are a generahzation of the Ehren- 
fest theorem of Hnear quantum mechanics to the nonhnear Schrodinger equation and particu- 
larized for the potential that concerns us [241 132| . 

This result has been discussed previously in many papers and is physically very interesting. 
It indicates that the evolution of the center of the solution is independent of the nonlinear 
effects and of the evolution of the rest of the moments and depends only on the potential 
parameters. 

4.3. Reduction of the general problem to the case I^^j — Vo,j — 0. We shall begin 
by stating the following lemma '33': 

Lemma 4.1. Let u{x,t) he a solution of (|4.1|1 with the initial datum u(x,Q) — uq{x). Then, 
the Junctions 



UR{x,t) =u(a:-i?(t),i)e'^(^'*) 



where 



(4.4a) 



e{x, t)^[x,R]+ J {R{t'),R{t')) - {R{t'),A{t')R{t')) 



dt' 



and 



dt^ 



AR = 



(4.4b) 



(4.4c) 



for any set of initial data i?(0),i?(0) G M" are also solutions of 14.1|l . 
Proof. All we have to do is substitute 14.4all , 



(|4.4b|l and (fOcIl in (gnj. □ 

One noteworthy conclusion is that given a solution of Eq. (|4.1|l we can translate it initially 
by a constant vector and obtain another solution. In the case of stationary states, defined as 
solutions of the form 



(4.5) 



which exist in the autonomous case (i.e. dX/dt — 0) and whose dynamics are trivial, this result 
implies that under displacements the only dynamics acquired are those of the movement of the 
center given by Eq. H4.4c(l . The coincidence of the evolution laws H4.3|l and H4.4c(l allows us to 
state the following theorem which is an immediate consequence of the above lemma. 

Theorem 4.2. Ifip{x,t) is a solution of 1)4. l|l with non-zero ifj or V^^, then there exists 

a unique solution u{x,t) — ip{x + {/^ ■(t)}j, t)e'^^^'*^ with 



e{x,t) 



E 



dt' 



such that = and Vq"^ = 

The important conclusion of this theorem is that it suffices to study solutions with Iij 
and Vqj equal to zero, as those that have one of these coefficients different from zero can be 
obtained from previous ones, by means of translation and multiplication by a linear phase in 
X. From a practical standpoint, what is most important is that Iij be null without any loss 
of generality, as then we can establish a direct link between the widths and the moments l2.j 
(see the discussion in the third paragraph of Section 
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4.4. Moment Equations. Assuming that all of the moments can be defined at any time 
t, we can calculate their evolution equations by means of direct differentiation. The results are 
gathered in the next theorem. 

Theorem 4.3. Let uq{x) be an initial datum such that the moments l2.j, Vij, Kj and J are 
well defined at t — 0. Then 



dl, 



2j 



dt 

dt 
dK, 



AK, 



dt 



D{p, t)dx, 



'^^lDip,t)A,dx^l 



dG{p,t) 
dt 



dx. 



(4.6a) 
(4.6b) 

(4.6c) 

(4.6d) 



where D{p,t) = G{p,t) - pg{p,t), p = \u{x,t)\'^. 

Proof. The demonstration of the validity of Eqs. H4.6|l can be carried out from direct 
calculations, performing integration by parts, and using the decay of u and Vm at infinity. 

To demonstrate ()4.6af) it is easier to work with the modulus-phase representation of u, 



u = p 



1/2 A4, 



(with p>0). Then 



dl 



2,i 



dt 



J x^jp = - j x]{Vp -Vq^ + pA-cj)) 

- J x]Vp -^(1) + j V{x]p) ■ V0 

~ j x']\7p-\7<i) + j x]\Ip-\7(t) + j \I{x])p-W(j) 



2 / Xjp 



50 
dxj 



We can also prove Ij4.6bp as follows 



d¥^ 
dt 



du 



dut 



Xj Ut 



du 



dxj dxj dxj 



i- . du i / 



\^jXj\u\ 



1 



= -2A 



/ X, 



2 

a;>p~2 

du 
dxj 



I— 

dxj 

du 

^^^^J^'d^, 

du dg 
dxj dxi 



Xlxl 1 u 



du 1 dAu 
^ dxj 2 dxj 
.n du 



I |2 ^9 

. _ du dAu _ dAu 
Au— u— h u 



dxj 



dxi 



dxi 



(4.7) 
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where c.c. indicates the complex conjugate. Operating on the above integrals we have 



du ^ - du dAu _ dAu \ 
Au- h Au- u— u— — 

oxj oxj oxj oxj J 



2 / Xi Vu ■ 



dVu 

dxj 



Vu ■ 



dVu\ 
dxj ) 



du 



dxi 



2 / |Vw|- 



du 



dxi 



and 



dg 

' dxi 



gp 



dG 

' dx-i 



gp 



G= D 



By substitution in (|4.7|l we arrive at the desired result. 



dt 



= -2A? / x]\u\'^ + 2 



du 



dxj 



'2Xp2,j + -2 I D. 



-2 D 



Let us now prove 



dKj 1 f d f du du 



2 J dt \ dxj dxj 



dut du du dut 
dxj dxj dxj dxj 



dxi 



2 2 1 

feXfe M - tgu 



\k=l 



du 
dxi 



c.c. 



then 

^^^ix" f X- (u—-u—] - - I ^(u—-u—\ 
dt 2 ^ J •' \ dxj dxj J 2 J dxj \ dxj dxj J 

i [ (d^d^ _ du_dAu\ _ _l^2^ f d / d<j>\ 

47 \ dxj dxj dxj dxj J 2 ^ J ^ dxj \ dxj J 

1 .2,r ^ [ ( d^4> ^dp dA 



and using the definition of G we obtain 

_ 1 ,2t. f d'^'t' f dG d(t> 
-dr-'2^=^'-^^ J ^^dx]^ J d^,d^, 
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Finally, to demonstrate Ij4.6d|l we proceed 



dJ 



dG{p,t) _ 
dt 
dG , 



dG_ 
dp 



dp dp . 



du 



Ut 



du 



Ut 



dG 
'dt 



_ , dG 
nut) + — 



-gV ■ (pV0) + — 



-VG- V0 

DA(t)- 



dG 
dt • 



g {uAu ■ 



A N 9G 
uAu) + — 
dt 



gpAcj) + — 



dG 
dt 



dG' 

~gVp ■ V(/. - .gpA0 + — 
= /(G-.p)A0 + /f 



A direct consequence of the theorem is 

Corollary 4.4. Letu{x,t) be a stationary solution of (|4.1|l . Then, 



i J D{p)dx. 



(4., 



5. Solvable cases of the method of moments. In this section we will study several 
particular situations of practical relevance in which the method of moments thoroughly provides 
exact results. 

5.1. The linear case g{p,t) — 0. In this case, Eqs. (j3.1d|l and l|4.6d|l tells us that 
J{t) — 0, for all t, and then the moment equations (|4.fi|l become 



dl: 



2j 



Vij, (5.1a) 



dt 
dt 



4Xj-2A^^/2j, (5.1b) 
dK, 1 , , ^ 

That is, in the linear case the equations for the moments along each direction ] of the physical 
space M" are uncoupled. This property was known in the context of optics for n — 2 and 
constant j84j . Here we see that this property holds for any number of spatial dimensions, 
time dependence A(i) and even for nonsymmetric initial data. 

5.2. Condition of closure of the moment equations in the general case. Equations 
(14. 6() do not form a closed set and therefore to obtain, in general, their evolution we would 
need to continue obtaining moments of a higher order, which would provide us with an infinite 
hierarchy of differential equations. Given the similarity among the terms that involve second 
derivatives of the phase of the solution in Eqs. H4.6(l . it is natural to wonder whether it would 
be possible to somehow close the system and thus obtain information about the solutions. 

From this point on, and for the rest of the section, we will limit ourselves to the case 
Aj(t) = A(t), J = 1, . . . ,n, which is the most realistic one, and which includes as a particular 
case the situation without external potentials Aj — 0. Let us define the following quantities 

n n n 

I - E ^2,,' V = 5] K^Y. (5-2) 

j=i i=i i=i 
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— = V, 

1 



2n 



Differentiating equations \b.2\ and using H4.6|l. we have 

dX 

H 
dV 
dt 

dJC _ 
dt ~ 
dJ _ 
H ~ 

If we add up equations H5.3c|l and 

dl 
'dt 
dV 
'dt 

d{JC + J) 



D{p, t)dx, 



Dip,t)A(b 

we arrive at 



D{p,t)A(l)dx, 

dG{p,t) 
dt 



dx. 



^ = V, 



= 4 



1 



dt 



-yv + 



D{p, t)dx 
dG 



dt 



dx. 



(5.3a) 
(5.3b) 
(5.3c) 
(5.3d) 

(5.4a) 
(5.4b) 
(5.4c) 



In order that equations 15.411 form a closed system, they must fulfill that — ^ Jg„ D{p,t)dx = 

J = /ign G(p, t)dx and that J-^„ ^^^^^dx can be expressed in terms of the other known quan- 
tities. The former condition requires that 



= 



-D{p,t) + G{p,t) 



dx ~ 



"\ ^/ ,\ " dG(p, t 



dx. 



As G{p, t) does not depend explicitly on a;, this condition is verified when 



dG{p,t) 
dp 



2 

1 + - 

n 



G(p,t) 



that is, if 



or, equivalently, if 



G(p,t)=5o(<)p'+'/", 



g(p,i) = .go(t)p2/", 



(5.5) 



(5.6) 



(5.7) 



where go{t) is an arbitrary function that indicates the temporal variation of the nonlinear term. 
Then 



/ 



dG{p,t) 
dt 



dx 



1 dgo 
90 dt 



G{p, t)dx 



1 dgg 
go dt 



J{t). 



(5.8) 



To close the equations it is necessary that .go(i) be constant in order to cancel the last term of 
this expression. Then, the nonlinearities for which it is possible to find closed results are 



9{p) ^ 9op""' = goH^'"-- 



(5.9) 



with gQ e M, remembering that in the case go < there may be problems of blowup in finite 
time. Fortunately, these nonlinearities for = 1, 2, 3 correspond to cases of practical interest. 
For instance, the case n = 1 with quintic nonlinearity has been studied in Refs. j35ll3fiir?7] and 
the case n — 2, with cubic nonlinearity, corresponds probably to the most relevant instance of 
NLS equation, i.e. the cubic one in two spatial dimensions |3(JI 131) . For n — 3 the nonlinearity 
given by Eq. H5.9|) appears in the context of the Hartree-Fock theory of atoms. 
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5.3. Simplification of the moment equations. Defining a new quantity £ — IC + J, 
Eqs. H5.4|l become 

^=V, (5.10a) 
^ =4£-2A^(i)X, (5.10b) 

f = -l/m (5.10c) 

These equations form a set of non-autonomous Uncar equations for the tliree averaged moments: 
£{t),V{t) and T{t). To continue our analysis, we note that 

Q = 2£I-^V^, (5.11) 

is a dynamical invariant of Eqs. (|5.10() . We finally define X = |Q| ^^'^J^/^, which is propor- 
tional to the mean width of u. A simple calculation allows us to corroborate that the equation 
that X{t) satisfies is 

Solving H5.12|l allows to calculate V and £ by simple substitution in H5.10|l . This equation 
is similar to that obtained in |30| for solutions of radial symmetry in the case n — 2 and 
g{u) = Here we find that it is possible to obtain a more general result for solutions without 
specific symmetry requirements, and for any combination of dimension and nonlinearity g{u) = 
\u\P satisfying the condition np = 4. The case with sgn(Q) ~ —1 corresponds to collapsing 
situations ITH|[3H1. In what follows we consider mostly the case Q > 0. 

Equation H5.12|l was studied by Ermakov in 1880 ^39; , although since then it has been 
rediscovered many times (see e.g. It is a particular case of the so-called Ermakov 

systems |41j ,42, ,43,, for which it is possible to give fairly complete results. Especially easy, 
though tedious to demonstrate is 

Theorem 5.1 (Ermakov, 1880). Let X{t) he the solution of (|5.12() with initial data X{Q) = 
Xo,X(0) = Xq. Then, «/xi(0 and X2{t) are solutions of the differential equation 

f^ + X\t)x^O (5.13a) 
satisfying the initial data xi(0) = Xq, Xi{0) = Xq and X2(0) = 0, X2(0) ^ then 



X{t)^\lxlit) + ^xlit) (5.13b) 



where w is the constant w — X1X2 ~ X2Xi- 

Equation (j5.13b|) is often called the principle of nonlinear superposition. Equation (|5.13a(l 
is the well known Hill's equation T4 which modelizes a parametrically forced oscillator, and 
which has been studied in depth. In the following, we shall study a couple of special situations 
in view of their physical interest. 

It is remarkable that the complex dynamics of a family of nonlinear partial differential 
equations can be understood in terms of a simple equation such as Hill's. 

If we suppose that the function A^(t) depends on a parameter e in the way A^(i) = 1 + Ae(t), 
Xe{t) being a periodic function with maximum value e (not necessarily small), there exists a 
complete theory that describes the intervals of values of e for which the solutions of Eq. (|5.13a|l 
are bounded (intervals of stability) and the intervals for which the solutions are unbounded 
(intervals of instability) '44' . 
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6. Applications. 

6.1. Dynamics of laser beams in GRIN media. When a laser beam propagates in a 
medium with a gradex refraction index (GRIN medium) with an specific profile quadratic in 
the transverse coordinates, the distribution of intensity u(x, y, z) in permanent regime is ruled 
by Eq. H4.1|l with g(p) = p and n = 2 (in the optical version of the equation i <-^- z), so that we 
are dealing with the critical case that we know how to solve. Although in principle it would be 
possible to design fibers with arbitrary profiles, technically the simplest way is to join fibers 
with different uniform indexes in each section. 

In this case, the modelization of the phenomenon is given by 

A^W41 '"["'"^"^ . (6.1) 

^' teiT,,Ta+n=T] ^ ^ 

Equation H5.13a|l with X{t) given by (|6.1() is known as the Meissner equation whose solution 
is trivial, given in each segment by a combination of trigonometric functions. 

The solutions to the Meissner equation can be bounded (periodical or quasiperiodical) or 
unbounded (resonant oscillations). In figure RTTl the two types of solutions are shown for a 
particular choice of parameters. 

As far as the regions of stability in the space of parameters are concerned, they can be ob- 
tained by studying the discriminant of equation H5.13a|l . defined as the trace of the monodromy 
matrix, that is 

Dia,b,Ta,n):^MT) + ^'2iT)- (6.2) 

where (j>i,<j>2 are the solutions of H5.13a|l satisfying the initial data (/)i(0) — l,<j)[{0) = and 
4>2{0) = 0, (/)2(0) — 1 respectively. 
In our case it is easy to arrive at 

(f>i (T) = cos{aTa) cos{bTb) - ^ sen(aT„) sen(5Tb), (6.3a) 

0^(r) = cos{aTa) coa{bTb) - - sen(aTa) sen(6Tb). (6.3b) 

a 
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Finally, the form of the discriminant is 

D{a, b, Ta, n) - 2 cos(ar, + bn) - ^° ~ sen(ar,) sen(6T6) (6.4) 

ab 

The Floquet theory for linear equations with periodical coefficients connects the stability of the 
solutions of (|5.13a(l with the value of the discriminant. The regions of resonance correspond to 
values of the parameters for which \D\ > 2, whereas if \D\ < 2 the solutions are bounded |44j . 
The equations D{a,b,Ta,Tii) = 2 and D{a,b,Ta,Ti,) = —2 are the manifolds that limit the 
regions of stability in the 4-dimcnsional space of parameters. In reality, defining a = aT, (3 = 
bT, Ta — ^T, Tb — (1 — 7)T the number of parameters is reduced to three 

D{j, a, /3) = 2 cos(a7 + /3(1 - 7)) - ~ sen{a-/) sen(/3(l - 7)), (6.5) 

ap 

Therefore, the isosurfaces D{'~f,a,(3) = 2 and _D(7, a,/3) = —2 can be visualized in three 
dimensions as is shown in figure 1^7^ 

The general study of the regions that appear in figure 1^?^ is complex, which leads us to focus 
below on a few particular cases. 

For example, in the case in which the two sections have the same length Ta — T^, the 
discriminant depends only on aT, bT and it is 

, f(a + b)T\ (aT-bTf faT\ f bT\ , , 

so that now the condition |Z?| = 2 determines curves such as those of figure I^TTT b) . 

The structure of the regions of resonance can be explored in more detail fixing the relative 
values of the coefficients, for example taking b = 2a. 

/3aT\ 1 ( aT 



D{a, T) = 2 cos j - - sen j sen(aT) (6.7) 

Defining the variable y — aT/2, the discriminant is a function D{y) [Fig. 16. 3j . The so-called 
characteristic curves are hyperboles of the form 2j/j?' = aT being y^j"-* , y^'* respectively the 
solutions of the algebraic equations f±{y) = 2cos3y — i sen y sen 2j/ =F 2 = 0. It is easy to 
demonstrate that the regions of resonance are contained between two consecutive zeros of /+ 
or /_ that can be obtained using any elementary numerical method. If we draw a as a function 
of 1/T, the regions of resonance are those marked in colors in figure I^^T b). Obviously the 
image is repeated due to the periodicity 27r of D{y) and there are only four basic regions 
of resonance (together with its harmonics) contained in the intervals (roots of /+ and /_): 
y G [0.84, 1.23] U [1.91, 2.3] U [3.98, 4.37] U [5.05, 5.44] (see figure Ela)). 

Another case of possible interest is that in which one of the fibers is not of GRIN type, that 
is, 6 = 0. Then the discriminant is given by the limit of equation H6.4|l when b ^ 0. 

L»(a,T) = -aT sen +2 cos (6.8) 

As in the previous case the only relevant parameter is y — aT/2, the regions of resonance on 
the plane a — T are hyperboles and the relevant quantities are the zeroes of /+ and /_ , which 
are given by 



f±{x) = —yseny + 2 cosy =F 2 = 0. 



(6.9) 
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Fig. 6.2. (a) Isosurfaces corresponding to D = —2 (in this range of values D < 2) for a limited range of 
parameters. The regions between the gray surface and the planes limiting the drawing are regions of resonance. 
Sections are shown for two particular values of fS where the bluish tones correspond to the regions of resonance. 
The color bar indicates the color corresponding to each level of D(7, a, 0) and the arrow indicates the color 
assigned to the isosurface D = —2. (b) The same as (a) but for a larger range of parameters. Isosurfaces 
D = 2 and D = —2 are shown in brown and green respectively. A section is shown for a particular value of j3 
with bluish and reddish tones corresponding to regions of resonances with D < —2 and D > 2, respectively. In 
this case the two values D = 2 and D = —2 are indicated by arrows on the color bar. 
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l/T 

Fig. 6.3. (a) First four regions of resonance in y (in color) for Ta = T/2,b = 2a, as a function of 
y = aT/2. (b) First five regions of resonance in the plane a — l/T. 

Now there is no exact periodicity in the positions of the zeros any more, but at least it is 
possible to estimate the location of those of high order. To do this we must bear in mind that 
for y big enough the dominating term in both cases is f±{y) — —yseny, so that the zeros will 
be given by y — nvr. It can be seen with a perturbative argument that the convergence ratio 
is of the order of 0{l/n). Writing y'^^ = nn + e^-* and substituting it in equations (|6.9|) it is 
found that 

4"' ~ (-i)"+ii±^. (6.10) 

This type of analysis can be extended to any restricted set of parameters. 

6.2. Dynamics of Bose-Einstein condensates. In the case of Bose-Einstein conden- 
sates, there has recently been great interest in the study of the dynamics of these systems in 
a parametrically oscillating potential. Recent experiments (see e. g. 0^1 E|) have motivated 
a series of qualitative theoretical analyses (the pioneer works on this subject can be seen in 
|48[ I49L ISUj although there is a great deal of posterior literature). 
In the models to which we refer, the trap is modified harmonically in time, that is 

X'^ (t) = l + ecos ujt (6.11) 

with e > — 1. Equation (|5.13a|l with A(t) given by (|6.11|) is called the Mathieu equation. For 
this equation it is possible, as in the case of the Meissner equation, to carry out the study of the 
regions of the space of parameters in which resonances occur. In the first place, for any fixed 
e, there exist two successions {w„}, {uj'n} with ujn,uj'^ — > such that if we take uj G (a;„, w^), 
equation (|5.13a(l possesses a resonance. In the second place, for fixed iv, the resonances appear 
when e is large enough. The boundaries of those regions are the so-called characteristic curves 
that cannot be obtained explicitly but whose existence can be demonstrated analytically, as in 
the previous section, by using the discriminant. In the case of the Mathieu equation, it can be 
proven that the regions of instability begin in frequencies uj ~ 2,1, 1/2, . . . , . . . |44j . 

As in the previous case, the resonant behavior depends only on the parameters, and not 
on the initial data. With respect to stability, the Massera theorem implies that if {e, uj) is 
in a region of stability, then there exists a periodic solution of (|5.13a|l . and by the nonlinear 
superposition principle such a solution is stable in the sense of Liapunov. 
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7. Approximate methods (I): Quadratic phase approximation (QPA). 

7.1. Introduction and justification of the QPA. Up to now, the results we have 
shown for the evolution of the solution moments are exact and in some sense rigorous. Unfor- 
tunatelly, in many situations of practical interest it is not possible to obtain closed evolution 
equations for the moments. In this section we will deal with an approximate method which is 
based on the method of moments. 

The idea of this method is to approximate the phase of the solution u by a quadratic function 
of the coordinates, that is 



where U{x,t) is a real function. 

Why use a quadratic phase? Although there is not a formal justification and we do not know 
of any rigorous error bounds for the method to be presented here, there are several reasons 
which can heuristically support the use of this ansatz for the phase for situations where there 
are no essential shape changes of the solutions during the evolution. First of all, when Eq. (|4.1|l 
has self-similar solutions, they have exactly a quadratic phase 021 ■ Secondly, the dynamics 
of the phase close to stationary solutions of the classical cubic NLS equation in two spatial 
dimensions (critical case) is known to be approximated by quadratic phases |E| Finally, 
to capture the dynamics of the phase of solutions close to the stationary ones, which have a 
constant phase, by means of a polynomial fit, the terms of lowest order are quadratic since the 
linear terms in the phase may be eliminated by using Theorem 14.21 

We would like to remark that a frequently used approach in the physical sciences to study 
the evolution of solutions of Eq. 1)1. l|l is the so called collective coordinate method (also called 
averaged lagrangian method or variational method) where an specific ansatz for the solution 
is proposed depending on a finite number (5) of free parameters u(x,t) — Lp{x,pi{t), ...ps{t)). 
The equations for the evolution of the parameters are then sought from variational arguments 
|2nil^. For NLS equations all commonly used ansatzs have a quadratic phase. Our systematic 
method provides a more general framework in which other methods can be systematized and 
understood. 

As we will see in what follows, the choice (|7.1|) allows to obtain explicit evolution equations 
and solves the problem of calculating the integrals of the phase derivatives in the equations 



7.2. Modulated pow^er-type nonlinear terms. Under the QPA, for modulated power- 
type nonlinearities g{p,t) ~ <?o(^)p^^^i p G M for which J^^^ D{p)dx = — pJ/2, the moment 
equations H4.6|l are 




(7.1) 



dl2,j 



(7.2a) 



dt 



4Kj - 2\p2,j + pJ, 
-\><]V,,,+pp,J, 



(7.2b) 



dt 
dK, 



(7.2c) 



dt 




(7.2d) 
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To these equations we must add the identity Vij = 4/3j/2j-, which is directly obtained by 
calculating Vij. Or, expressed otherwise 



P. - (7.3) 



Let us now consider the simplest case of solutions with spherical symmetry with = A(t), 
j = 1, . . . , n, for which (f){xi, a;„) = (3{t) (x^ + ... + a;^). Using the same notation as in Eqs. 
1)5. 2|l the moment equations become 



dJ 



V, (7.4a) 



dV . f ^ 



dt 

dK 1 



4(aC + ^j) -2A2x, (7.4b) 
A^V + np/3J, (7.4c) 



dt 2 

^ = -np/3J+l%J (7.4d) 
dt go dt 

with V = 4/3X. Despite the complexity of the system of equations (|7.4(l it is possible to find 
two positive invariants 

Qi = 2/CX-VV4, (7.5a) 

Q2 = ^r^P'\j. (7.5b) 

The existence of these invariants provides J as a function of X, which allows us to arrive at an 
equation for X — X^l"^ 

'^ + mx^^ + go{t)^^. (7.6) 

Again we obtain a Hill's equation with a singular term. Note that in the case n — 3,p — 2 
we have a quartic term in the denominator, which corresponds with the type of powers that 
appear in the equations which are obtained in the framework of averaged Lagrangian methods 

1201. 

The quadratic phase method provides reasonably precise results that at least desribe the 
qualitative behavior of the solutions of the partial differential equation. Using several numerical 
methods, we have carried out different tests especially in the most realistic case = 6 in 
(|4.1|) . For example, in figure ITTI the results of a simulation of H4.1|l with n ^ 3, p — 2, 
A'^(t) = 1 + 0.1sin(2.8i) and ga — 10 are presented for an initial datum uo{x) = e^^ ^^/n^^^. 
In this case the simplified equation H7.6|) predicts quasi-periodic solutions which is what we 
obtain when resolving the complete problem. 

In Fig. 17. 21 we show the results for A^(t) = 1 + 0.1 sin(2.1t) for which H7.6|l predicts resonant 
solutions. Again, the results of the two models are in good agreement. 

Another interesting application of the quadratic phase approximation method is the case of 
cubic nonlincarity, g{p,t) = go{t)p — 9o(t)\u\'^, without potential X{t) = 0. In that situation 
Eq. (|7.6|l becomes 

d^X Qi Q2 
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Fig. 7.2. Results of the simulation of 14.11 in three dimensions on a grid of 64 X 64 X 64 with At = 0.005 
for p = 2, A^(i) = 1 + 0.1 sin(2.1t). X^{t) obtained numerically from the solutions on the 3D grid is shown. 



where the conserved quantities are 

Qi = 2/CX-VV4, (7.8a) 



n 



Q2 - — (7.8b) 

90 



This model describes the propagation of Ught in nonhnear Kerr media as weh as the dynamics 
of trapless Bose-Einstein condensates. In this situation the previous equations are used to study 
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the possibility of stabilizing unstable solutions of the NLS equation by means of an appropriate 
temporal modulation of the nonlinear term, that is, by choosing a suitable function go{t) thus 
providing an alternative to more heuristic treatments jSlL 1521 . More details can be seen in 
Ref. ISl. 

7.3. Closure of the equations in other cases. We have just seen that the quadratic 
phase approximation method allows us to close the moment equations in the case of power-type 
nonlinear terms. Following those ideas we have managed to close the equations in more general 
cases. 

We start from the evolution equations for the mean moments (|5.4|l . that after performing 
the quadratic phase approximation become 

f^V. («a, 

^ =4/C-2A2X-2n / D{p,t), (7.9b) 
at Jjjn 

^ = -^A^V - 2n(3 j ^ D{p, t), (7.9c) 

'^-2np[ DM^f (7.9d) 
at Jk" Jr" ot 

and V = 4/31. 

The idea to close the previous equations is to calculate the evolution of Jg„ D{p, t)dx, which 
is the term that prevent us from closing the equations, and try to write this evolution in terms 
of the moments. Let us define a new moment as 

T{t)= I D{p,t)dx = J- [ p^^^dx. (7.10) 
Jr" Jr" op 



Then, the evolution equations are 



4:JC-2X^I ~2nT, (7.11b) 



§ = V, (7.11a, 

dV 

m 

^ = -\\^V-2nPT, (7.11c) 

^ = 2nPT+ [ ^dx, (7.11d) 
dt 7r„ dt 

together with the evolution of T 

dJ" „ . „ . /■ 2d^G , f dG , f d dG 



,^ 2n(3F + 2nP P^dx + ^dx ~ P^.ir'^x. (7.12) 
dt 7r" dp^ 7r" dt 7r„ dt dp 

To try to close the system of equations (|7.11|) - H7.12|I we impose that /jj„ p^^^dx is a linear 
combination of and J 

f d^G f f dG 

/ p^-^r-irdx = arJ^ + ajj = {ajr + aj) I Gdx - ayr p——(lx, (7.13) 
Jr" op Jk" Jm" cip 

where ajr and aj are two arbitrary constants. Then, G must verify 

.a^G dG 



0. 
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Therefore, if the nonUnear term g{p) in the NLS equation is such that G{p) verifies Euler's 
equation 

the evolution equations wiU close. In that case we can write G{p, t) — ga{t)Gi{p), where 5o(i) 
is an arbitrary function which indicates the temporal variation of the nonlinear term and Gi {p) 
satisfies Eq. (|7.14(l . So 

dt dt 7r,. go at 

d dG dgo f dGi 1 dgo . , . 
dt dp dt Jr,. dp go dt 

and the moment equations are written as 

f.V, (7^16.) 

dV 
'dt 

^^~\\^V~2nPT, (7.16c) 

2nPT+-'^J, (7.16d) 



AK,-2YX-2nT, (7.16b) 



dt go dt 
^ = 2n/3(l + a^)^ _ , 



'^"^ 2n/3(l + a^)JP + 2n/3ajJ+ — ^JP. (7.16e) 



By solving Eq. (|7.14|l we obtain specific nonlinear terms for which the quadratic phase ap- 
proximation allows us to write closed equations for the moments. Depending on the parameter 
(5 = (1 + ajr)^ + 4a J there exist three families of solutions 

'CipP+ +C2pP-, (S > 0, 
Gi{p) = \ Cip« + C2/9«logp, 5^Q, (7.17) 
Ci/o^ cos(/logp) + C2p^sin(/logp), (5 < 0, 

where p± = ((1 - a^) ± 5^'^) /2, i? = (1 - a^) /2, / = \5\^ /2. 

The most interesting case for applications is (5 > the nonlinear term being of the form 

gi(p) = fcipP+-i + fc2P^'-\ (7.18) 

where fci and ^2 are arbitrary constants and p+ and p- are defined through the relations 

a^F = 1 - - (7.19a) 
aj = -(p+-l)(p--l). (7.19b) 

Eq. H7.18|l implies that the quadratic phase approximation allows to close the moment equations 
for nonlinear terms which can be written as a linear combination of two arbitrary powers of 
\u\. 
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As in the previous subsection it is possible to find some invariant quantities, namely 

gi = 2/CX-VV4, (7.20a) 
, n 1"+ 

/+ 90 
, n 

/+ 50 

where 



Q+ = C- (J + /+^), (7.20b) 

C- (J + /-^), (7.20c) 



These conserved quantities allow us to write a differential equation for the dynamical width 

X(i) = jl/2 



_ A^W^-^+^^o(t)^^2^-^5^j. (7.22) 

The most interesting kind of nonlinearity in the form of Eq. (|7.18|) is the so-called cubic- 
quintic nonlinearity for which .go(^) = Ij 5i(p) = ^iP + ^29^ = fcil^P + A:2|m|^. Then, we have 
p+ - 1 = 2, p_ - 1 = 1, = -4, aj - -2, a+ - 1, /+ = 1, a_ - 1/2, /_ - 1/2, C = 2. The 
invariant quantities are 

Qi = 2/CX-VV4, (7.23a) 
Q+ = 2nX"(J + J^), (7.23b) 
Q_ = 2nX"/2( J + JC-/2), (7.23c) 

and the equation for the width is 

^ + A (t)X-— + (7.24) 

These equations contain a finite-dimensional description of the dynamics of localized solutions 
of the model and are similar to those found under specific assumptions for the profile u{x,t) 
(see e.g. [55. .56,. 57,). The main difference is that the method of moments allows to obtain the 
equations under minimal assumptions on the phase of the solutions and that depend on general 
integral quantities related to the initial data Qi, Q+, Q„. This is an essential advantage over 
the averaged Lagrangian methods used in the literature for which the specific shape of the 
solution must be chosen a priori (see also |2()l I45| ). 

8. Approximate methods (II): Thomas-Fermi Limit. 

8.1. Concept. In the framework of the application of the NLS equations to Bose-Einstein 
condensation problems (thus for nonlinearities of the form g(p) = gop), the Thomas-Fermi limit 
corresponds to the case go 3> 1 (note that this is only one of the many different meanings of 
"Thomas- Fermi" limit in physics). 

Usually, what is pursued in this context is to characterize the ground state, defined as the 
stationary solution of the NLS equation given by Eq. H4.5|l with fixed norm having minimal 
energy E. It is also interesting to find the dynamics of the solutions under small perturbations 
of the ground state solution. 
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8.2. Physical treatment. Let us consider the problem of cliaracterizing tlie ground state 
of H4.1|) . The usual "physical" way of dealing with this problem consists of assuming that if the 
nonlinear term is very large then it would be possible to neglect the Laplacian term in (|4.1|l 
(!) and to obtain the ground state solution as 



(fiTpix) 



\ 



90 




(8.1) 



The value of fi is obtained from the condition of normalization ||iy9Ti?||2 = 1. This procedure 
provides a solution without nodes which is then argued to be an approximation to the ground 
state. 

This method is used in many applied works, but unfortunately it is not even self-consistent. 
Near the zero of the radicand of 1)8. l|l the approximation obtained has divergent derivatives, 
which contradicts the initial hypothesis of "smallness" of the Laplacian term. Although several 
numerical results can be obtained using this approximation, its foundation is very weak. 

In order to understand the problem better, we rewrite 1)4. l|l making the change of variables 
K = n/go, rj = x/^/g^, ip{r]) = ip (x/^/g^), to give us the equation 



with e = 1/.90- It is evident that e^Aip is a singular perturbation whose effect may not be 
trivial. 

8.3. The method of moments and the Thomas-Fermi Umit. What can be said for 
the case of power-type nonlinearities in the limit 5 3> 1 on the basis of the method of moments? 
Before making any approximations we write an evolution equation for T as follows. For the 
sake of simplicity, though it is not strictly necessary, we will consider the case of Xj = A for 
j — 1, n and study the equations for the mean values (|5.2|) . 
First, we write Eqs. H5.3a|l and (|5.3b)l 

dl 

- = V, (8.3a) 



4:[IC + ^J] -2X^1 = {4:- 7ip)K. + npn- [2 + ^) X-'I. (8.3b) 



where Ti is the conserved energy. Combining H8.3a|l and (|8.3bp we arrive at 

A^X = (4 - np)IC + npH. (8.4) 



Equation H8.4|) is exact. 

The fact that the energy functional E reaches a miminum over ipo implies, by Lyapunov 
stability, that initial data uo{x) — ipo{x)+eS{x) close to the ground state must remain proximal 
for sufficiently small values of e. 

The only approximation needed to complete our analysis is to assume that when <? ^ 1 then 
J ^ K. for the ground state. Notice that this is a much more reasonable assumption than 
the direct elimination of the second derivative in the evolution equation. Thus, the energy 
conservation and the previous considerations allow us to affirm that J{t) ^ )C(t) for all times. 
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Although these facts can be used to write exphcit bounds for /C, as a first approximation 
and just in order to show the power of these ideas we can simply take /C ~ 0. Under this 
approximation we have 

g+(2+f)A^I«.npH. (8.5) 
whose solutions can be obtained explicitly as 

The equilibrium point of Eq. H8.4|) (corresponding to ^ = i? = 0) gives us the "size" of the 
ground state as a function of the physical parameters. Also the frequency of the oscillations 
around the equilibrium point is immediately obtained from Eq. (|8.6|l 



n^X^2+f. (8.7) 

We have performed numerical simulations of the partial differential equations H4.1|l to verify 
this prediction. Specifically, taking g = 5000, 20000, A = 1 and initial data of the form 
Uo{x) = (foii^ + e)x)/y/l + e for e = 0.01 and e = 0.02, we find a numerical frequency of 
f^num = 2.26 which is in excellent agreement with the value provided by our Thomas- Fermi 
formula f^xF = \/8 = 2.24. 

9. Summary and Conclusions. In this paper we have developed the method of mo- 
ments for Nonlinear Schodinger equations. First we have found the general expressions of the 
method and classified the nonlinearities for which it allows a closed explicit solution of the 
evolution of the moments. We have also discussed several applications of the method such as 
the dynamics of Kerr beams in nonlinear stratified media and the dynamics of parametrically 
forced Bose-Einstein condensates. 

Approximate techniques based on the method of moments are also discussed in this pa- 
per. In particular, the quadratic phase approximation is also developed here and applied to 
different problems, such as the writing of simple equations describing the stabilization of soli- 
tonic structures by control of the nonlinear terms and the dynamics of localized structures 
in cubic-quintic media. Finally, we have also studied the moment equations in the so-called 
Thomas- Fermi limit. 
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